Growth, microstructure, and failure of crazes in glassy polymers 



Jorg RottlejQ and Mark 0. Robbins 
Department of Physics and Astronomy, The Johns Hopkins University, 
3400 N. Charles Street, Baltimore, Maryland 21218 
(Dated: February 2, 2008) 

We report on an extensive study of craze formation in glassy polymers. Molecular dynamics 
simulations of a coarse-grained bead-spring model were employed to investigate the molecular level 
processes during craze nucleation, widening, and breakdown for a wide range of temperature, poly- 
mer chain length N, entanglement length N e and strength of adhesive interactions between polymer 
chains. Craze widening proceeds via a fibril-drawing process at constant drawing stress. The exten- 
sion ratio is determined by the entanglement length, and the characteristic length of stretched chain 
segments in the polymer craze is JV e /3. In the craze, tension is mostly carried by the covalent back- 
bone bonds, and the force distribution develops an exponential tail at large tensile forces. The failure 
mode of crazes changes from disentanglement to scission for N/N e ~ 10, and breakdown through 
scission is governed by large stress fluctuations. The simulations also reveal inconsistencies with 
previous theoretical models of craze widening that were based on continuum level hydrodynamics. 

PACS numbers: PACS: 81.05.Lg, 62.20.Fe, 83.10.Rs 



I. INTRODUCTION 

The failure of glassy polymers such as polystyrene 
(PS) or polymethylmethacrylate (PMMA) under exter- 
nal stresses occurs either through shear deformation or 
through crazing 0, 0] ■ While shear yielding occurs essen- 
tially at constant volume, crazing has a strong dilational 
component, and the volume of the material increases to 
several times its original value before catastrophic frac- 
ture occurs. Crazing is a failure mechanism unique to en- 
tangled polymeric materials and usually precedes a crack 
tip (see Fig.0. The fundamental and technological im- 
portance of crazes is that they are in part responsible for 
the large fracture energy G c of polymer glasses El 
that makes them useful load-bearing materials. They 
control the crack tip advance and require a large amount 
of energy dissipation up to the point of catastrophic fail- 
ure. Crazes can reach several fj,m in width and consist of 
an intriguing network of fibrils and voids that spans the 
entire deformed region. 

Despite the frequent appearance of crazes, there is still 
comparatively little theoretical understanding about the 
conditions and mechanisms of craze nucleation, growth 
and ultimate breakdown P, 0, IE 0> ■ In this paper, 
we present an extensive set of nonequilibrium molecular 
dynamics (MD) simulations that address these various 
phenomena. In this approach, polymers are modeled on 
a coarse-grained scale that takes into account van-der- 
Waals (vdW) and covalent interactions without specific 
reference to chemical detail. The effect of chain length, 
temperature T, widening velocity v and vdW interaction 
strength on the craze structure can be studied over a 
wide range of parameters. The molecular simulations 
allow insight into microscopic details not accessible to 



'Electronic address: rottler@pha.jhu.edu 



craze widening 











Crack 


d ~ 1jjm 


Craze I 


B craze 


^^^^^^ tip advance 
dense polymer glass 



craze length I ~ 0.1 mm 



FIG. 1: Craze fracture of glassy polymers. The craze is a 
deformed region (shaded) that grows in width and length un- 
der an applied vertical stress S. Its density is reduced with 
respect to the undeformed polymer by a constant extension 
ratio A. Characteristic values for width d and length I are 
indicated. During growth, S acts perpendicular to the sharp 
interface between undeformed polymer and craze. Also shown 
is an advancing crack tip from the left that breaks the craze. 
Small representative volumes of each region are studied with 
molecular simulations. 



experiments and offer an opportunity to test and develop 
theoretical models of crazing. 

A fundamental limitation on molecular level treat- 
ments is of course the finite system size. The largest 
volumes accessible at present are ~ lOOnm 3 , while the 
craze spans many fim. We are thus limited to a study of 
craze widening in a small representative region and can- 
not include e.g. the entire crack tip. Craze tip advance 
processes are beyond the scope of the present work. 

Several aspects of craze physics have already been 
adressed with simulations in previous papers. Baljon 
and Robbins 0,13 demonstrated the importance of chain 
length for the onset of craze growth. Rottler et al. stud- 
ied the elastic properties and fracture stresses of fully 
evolved crazes and used them in combination with linear 
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fracture mechanics to calculate the macroscopic fracture 
energy of glassy polymers that fail by crazing 6] . Rot- 
tier and Robbins also investigated how polymer entan- 
glements affect the craze structure on a microscopic level 
and argued that they "jam" the expansion of the glass 
under tension jlfij . 

This paper extends the previous work and is organized 
as follows. In Section [n] we briefly summarize the key 
experimental observations and review existing theoret- 
ical models of crazing. Section IIIII gives the technical 
details of the molecular models used in this study. We 
then analyze results for craze nucleation (Section II Vt . 
growth (Section EJ), microstructure (Section IVI|) . and 
failure fSection IVII|) . and compare our findings to pre- 
vious models and experiments. Final conclusions are of- 
fered in Section IVIIII 



II. CRAZE PHENOMENOLOGY AND THEORY 
A. Experiments 



zone h (see also Fig- EJ) is usually between (D) and (Do). 
Craze widening is a steady-state process, in which a con- 
stant "drawing stress" S ranging between 20 and lOOMPa 
is applied. Typical experimental values are 35 MPa 
(polystyrene) [lj] and 70 MPa (polymethymethacrylate) 
|38j . The value of S is of the same order as the shear 
yield stress of the polymer, and is found to increase with 
the entanglement density. 



B. Theory 

A theory of crazing has to explain the molecular ori- 
gin of the craze structure and the interdependencies of 
the various quantities measured in experiments. Despite 
a wealth of experimental data on crazing, there is cur- 
rently no theoretical description that addresses all as- 
pects of craze physics. The following models have been 
proposed to explain the extension ratio A and the rela- 
tionship between fibril spacing (Do) and drawing stress 
S. 



Crazes have been studied experimentally for more than 
30 years [H IH IH d The techniques most com- 
monly used to analyse the craze structure are transmis- 
sion electron microscopy (TEM) (l7lll8|| . low angle elec- 
tron diffraction (LAED) and small angle x-ray scattering 
(SAXS) [ijj UlA H3 ■ Comprehesive reviews of theoretical 
and experimental results have been presented by Kramer 
and Berger HHHl] and Creton et al. 

The density in the undeformed polymer pi and craze 
Pf is obtained from TEM measurements. The increase 
in volume during craze formation, or extension ratio A = 
Pi/pf, is found to have a characteristic value for a given 
polymer that is independent of molecular weight. Typical 
values of A for different polymers range from two to seven. 

Real space images of the craze show that the polymers 
are bundled into fibrils that merge and split to form an 
intricate network. The fibrils are highly aligned with the 
applied tensile stress and vary in diameter and length. 
However, this complex structure is normally idealized as 
a set of uniform vertical cylinders connected by short 
cross-tie fibrils The characteristic fibril diameter 

(D) and separation (Do) (see Fig- EJ are then determined 
from a Porod analysis of scattering experiments (see Sec- 
tion 01 . Measured values r ang e between 3 — 30 nm for 
(D) and 20-50 nm for (D ) [Illllllllll For example, 
for polystyrene one obtains (D) ~ 6nm and (Do) ~ 20nm 

Nucleation of crazes |l2( occurs preferentially near de- 
fects in the polymer. These produce large local tensile 
stresses that lead to the formation of microvoids that 
evolve into a craze. Once nucleated, the craze grows 
in length and width (see Fig. . It is well established 
|15| that the craze widens by drawing material from the 
dense polymer into new fibrils. This deformation is con- 
fined to a narrow active zone at the interface between 
the dense polymer and craze. The width of the active 



1. The extension ratio A 

The extension ratio A has been successfully explained 
by a simple scaling argument, that relates A to the micro- 
scopic entanglement network in the polymer glass. En- 
tanglements arise in dense polymeric systems from the 
topological constraints that the chains impose upon each 
other. The mobility of the chains is greatly restricted, be- 
cause they cannot pass through each other. The starting 
point for the present argument is the assumption that the 
glass inherits these entanglements from the melt, where 
an entanglement molecular weight is given by the plateau 
modulus under shear, . 

M e = plRT/SG%\ (1) 

This result can be derived from the microscopic tube 
model [24(, which relates the rheological response of the 
polymer melt to the deformation of a tube to which 
the polymer chain is confined. With repeat units of 
weight Mo, one can define a typical number of steps 
N e = M e /Mo (entanglement length) between entangle- 
ments along the polymer backbone. 

These entanglements are assumed to act like perma- 
nent chemical crosslinks during crazing, which implies 
that the expansion ends when segments of length iV e 
are fully stretched. The initial separation of entangle- 
ment points is di — (IploNe) 1 / 2 , according to standard 
random walk (RW) scaling, where lo is an elementary 
step length and l p the persistence length. The length 
of this segment rises from di to a maximum final length 
df = A max di = NJ Q , and thus 

A max = (NJo/lp) 1 ' 2 - (2) 
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Experimentally, Eq. J5J is well confirmed, but Section 
IV PI shows that the picture motivating this expression is 
oversimplified. 



2. The drawing stress S 

The value of the drawing stress S has traditionally 
been related to the craze microstructure ((D), (Do)) via 
capillary models In these models, the polymer in 

the active zone is treated as a viscous fluid with a surface 
tension T and a viscosity 77. Figure shows an idealized 
picture of the craze geometry, where craze formation is 
modelled as the propagation of void fingers with a char- 
acteristic spacing (Do) into the strain-softened fluid. The 
applied stress S required to advance the interface has a 
dissipative contribution arising from a suitable flow law 
(e.g. power-law fluid) and an energy penalty contribu- 
tion due to the surface tension. The tension is S in the 
polymer glass and the Laplace pressure 2T/((Dq)/2) at 
the ceiling of the finger, where (Dq)/2 is the character- 
istic radius of curvature (see Fig. By estimating the 
width of the active zone as h ~ (Do)/ 2, Kramer calcu- 
lated a stress gradient between glassy polymer and the 
finger void ceiling, 
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Since Vcr is proportional to the interface velocity, he then 
predicted that the system will select a value of 



(D ) ~ 8T/S, 



(4) 



which maximizes the stress gradient between finger ceil- 
ing and bulk polymer and thus will lead to the fastest 
propagation velocity of the fingers. 

More recently, Krupenkin and Fredrickson |7J have for- 
mulated a theory of craze widening that is similar in spirit 
to Kramer's arguments and also equates the craze widen- 
ing stress with a viscous and a surface tension contribu- 
tion. However, these authors suggest a different interpre- 
tation of r. They introduce an effective surface tension 
that begins to rise above the vdW value when the finger 
radius rises above the rms spacing between entanglement 
lengths di . This ansatz is motivated by the idea that ex- 
panding the random walk between entanglements gener- 
ates an additional energy penalty. An upper bound to 
r is provided by the energy required for chain breaking, 
which sets in once the finger radius exceeds the maxi- 
mum elongation between entanglement points, loN e . By 
minimizing the finger propagation stress, they conclude 
that the fibril spacing will always be 



Do 



(5) 



independent of surface tension. In their model, the fibril 
spacing is determined exclusively by the entanglement 
network. 




FIG. 2: Surface tension model of craze widening. Void fin- 
gers with characteristic spacing (-Do) propagate into a strain- 
softened layer of polymer fluid of width h, leaving behind 
fibrils of a characteristic diameter (D). The externally ap- 
plied stress S acts perpendicular to the fluid-glass interface. 
The characteristic radius of the finger caps is on the order of 
(D )/2. (See Refs. EHH for an analogous figure.) 



III. SIMULATIONS AND MOLECULAR 
MODELS 

We study craze formation by performing molecular dy- 
namics simulations of a standard coarse-grained polymer 
model [25| , where each linear polymer contains N spher- 
ical beads of mass m. Models of this kind have a long 
tradition in polymer research and have verified theories of 
polymer dynamics |24j in the melt. They have recently 
been employed by other researchers to study failure in 
network polymer adhesives |27| and end-grafted polymer 
chains between surfaces p^. 

In this bead-spring model, van der Waals interactions 
between beads separated by a distance r are modeled 
with a truncated Lennard- Jones potential: 

Vtj(r) = 4n [(a/r) 12 - (a/rf - (a/r c ) 12 + (a/r c f] 

for r < r c , where uo ~ 20 — 40 meV and a ~ 0.8 — 1.5 nm 
are characteristic energy and length scales 26] . A simple 
analytic potential [28| 



V hT (r) = -h(r - R f(r - RJ 



(7) 



is used for covalent bonds between adjacent beads along 
the chain. The form of this potential was chosen to al- 
low for covalent bond breaking, which is not possible 
with other standard bond potentials such as the popu- 
lar FENE potential [2f| . Bonds are permanently broken 
when r exceeds Rq = 1.5 a. The constant i?i = 0.7575a 
was chosen to set the equilibrium bond length lo = 0.96 a, 
which is the "canonical" value for the bead-spring model 
with the FENE potential |2j|. This allows us to use 
results from previous studies, most importantly the en- 
tanglement length. The constant k\ determines the ra- 
tio of the forces at which covalent and van der Waals 
bonds break. We find that this ratio is the only impor- 
tant parameter in the covalent potential and set it to 100 
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based on data for real polymers 27, 28], which implies 
k\ = 2351wo/a 4 . Tests with other analytical forms of 
the bond potential showed no appreciable impact on our 
results as long as the bonds break before the chains can 
pass through each other. 

In order to vary the entan glem ent length, we include 
a bond-bending potential [28ll29j] 



N-l 



Vs = b J2 1 



i=2 



(fi-i - rj) ■ (n - f i+ i) 
\(n-i - ri)\\(fi - f i+ i)| 



(8) 



that stiffens the chain locally and increases the radius of 
gyration. Here, n denotes the position of the ith bead 
along the chain, and b characterizes the stiffness. We will 
consider two cases here referred to as flexible (b = 0) and 
semiflexible (b = I.5ito) polymers. The corresponding 
entanglement lengths are » 70 and iV| fl 30 beads, 
respectively pllMl^. 

We consider three temperatures T = O.OIito/fcs, T = 
O.luo/fcs and T = 0.3uo/ks, where the last tempera- 
ture is close to the glass transition temperature. The 
amount of adhesive interaction between beads is varied 
by changing the range r c of the LJ potential from 1.5a 
to 2.2a. 

The equations of motion are solved using the veloc- 
ity Verlet algo rithm wi th a timestep of dt = 0.0075 tlj, 
where tlj = \J ma 2 Juq is the characteristic time given by 
the LJ energy and length scales. Periodic boundary con- 
ditions are employed in all directions to eliminate edge ef- 
fects. The temperature is controlled with a Nose-Hoover 
thermostat (thermostat rate Itjj 1 ), and the thermostat 
is only employed perpendicular to the direction of craze 
growth. Simulations with a Langevin thermostat showed 
no appreciable difference between the two methods. 

In all simulations of crazing, an initial isotropic state 
in a cubic simulation cell of edge length L is created using 
standard techniques 26] . Polymer chains are constructed 
as ideal RWs with a suitably chosen persistence length 
lp. l p is fixed by matching the radius of gyration of the 
chains to the equilibrium value in the melt, and the values 
are lp — 1. 65a and lp a = 2.7a for flexible and semiflex- 
ible chains, respectively. Subsequently, the interaction 
potentials are imposed and the system is cooled at con- 
stant volume from a melt temperature T m = I.Suq/Icb, 
to the desired run temperature. 

All runs begin at zero hydrostatic pressure. Strains eu 
are then imposed by rescaling the simulation box periods 
Li and all particle coordinates proportionately [3(j ■ This 
allows arbitrary stress states to be studied in Sectiou lIVI 



IV. CRITERIA FOR CAVITATION AND CRAZE 
NUCLEATION 



The loading conditions on the polymer glass determine 
whether it will fail initially by shear yielding or the for- 
mation of voids and cavities. In general, strong triaxial 




FIG. 3: Octahedral shear stress r^ ct at yield as a function of 
pressure p at two different temperatures T = 0.3uo/fcs (open 
symbols) and T = 0.01tto/&s (filled symbols). The solid lines 
are fits to Eq. © and the dashed lines show the onset of cav- 
itation. Values of a are indicated for the two temperatures. 
Also drawn is a dotted line through the transition points that 
separates the regions of shear and cavitational failure. Here 
yield is associated with the strain where r oc t peaks. 



tensile stresses will favor cavitation. Cavitation and craz- 
ing are closely related, because crazes usually require the 
initial formation of microvoids ^5|. We therefore first 
address the initial failure of the polymer glass through 
either shear yielding or cavitation, and later discuss the 
formation of crazes. 

The loading conditions that lead to shear yielding in 
many experimental polymers |3lLl32| are most accurately 
represented by the pressure-modified von Mises yield cri- 
terion. It is formulated in terms of simple stress invari- 
ants, the hydrostatic pressure p — — (a\ + a 2 + 03)/3 
and the deviatoric or octahedral shear stress r oct = 

((0-1 - cr 2 ) 2 + (02 - 0-3) 2 + (0-3 - cn) 2 ) 1/2 /3, where the 
<7j denote the three principal stress components. The 
pressure-modified von Mises criterion states that yield 
will occur at an octahedral yield stress T^ ct given by 



r + ap, 



(9) 



where To is the yield stress at zero hydrostatic pressure 
and a is a dimensionless constant. Its physical motiva- 
tion is that the elastic free energy stored in shear defor- 
mation is proportional to T 2 ct and failure should occur 
when this energy exceeds a threshold that rises slowly 
with p. 

In ref. I'S-il we examined a much larger range of stress 
states than in previous experimental studies and showed 
that the pressure-modified von Mises criterion provides a 
good description of shear yield in our bead-spring model. 
Data for two extremal temperatures are replotted in 
Fig.|3]along with solid lines showing fits to Eq. Shear- 
yield was observed to the right of the dot-dashed line, and 
these data points follow Eq. © quite accurately. To the 
left of the line cavitation was observed. The deviation 
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from the von Mises fits is very sharp, and r oc t quickly 
drops to zero. The values of r^ ct where cavitation occurs 
are well described by a straight lines 

T o c ct = r c + a c p (10) 

with new constants Tq and a c . This new "cavitation cri- 
terion" can be motivated in analogy to the von Mises 
criterion by assuming that the elastic free energy Fy as- 
sociated with volume changes must reach a critical value 
for cavitation to occur. Fy is proportional to p 2 , which 
gives a criterion of the form p — p . One can then assume 
that shear components in the stress tensor aid cavitation 
in a linear fashion, i.e. p = po + T oct /a c , which can be 
rearranged to give Eq. i|10|) with Tq = a c po. 

No clear experimental consensus exists about the stress 
state required for crazing, partly because of the impor- 
tance of surface defects in nucleating crazes. However, 
several criteria for craze nucleation were proposed almost 
30 years ago. They all try to take into account the crit- 
ical role of tensile stress components. Sternstein et al. 
[35| suggested a craze yield criterion of the form 

Tmax = gk* ~ fjlmax = A + B/p, (11) 

where A and B are constants that depend on tempera- 
ture. With respect to our criterion Eq. I |10| l. p has been 
replaced by 1/p and T oct by the largest difference between 
any two stress components. Bowden and Oxborough 1] 
formulated a similar criterion, where r max is replaced by 
(7i — VGi — v<73 and v is Poisson's ratio for the polymer 
glass. This expression is another possibility to describe 
the shear components of the stress state, and it reduces 
to T max when v — 1/2 and 02 — 03. The Sternstein 
and Bowden and Oxborough expressions could in prin- 
ciple also be fitted to the rather narrow range of pres- 
sure in Fig. where cavitation occurs. However, we are 
unaware of a convincing physical motivation for the 1/p 
term, which leads to obvious analytical problems at small 
p. In addition, the experimental results that motivated 
Eq. ljTT|) are sensitive to surface defects [T^ |. 

The above considerations pertain to the initial mode 
of failure of the polymer glass at strains typically less 
than 10%. However, crazing is a large strain deforma- 
tion with strains of several hundred percent. Although 
we find voiding to be a neccessary precursor to crazing, 
it is not guaranteed that a loading state that leads to 
cavitational failure according to Eq. Q10[) will ultimately 
produce stable crazes. Likewise, we have observed that 
an initial failure through shear deformation can still lead 
to later void formation and crazing. One should thus 
strictly call Eq. l|10fl a cavitation failure criterion and 
not a craze yielding criterion. 

V. GROWTH OF CRAZES 

In order to induce crazing, we enforce cavitation by ex- 
panding the periodic simulation box in the z-direction at 



constant velocity while maintaining the simulation box 
periods in the perpendicular x — y plane. This leads to 
an initial stress state where all three principal stresses 
are tensile. The initial voids formed during cavitation 
expand upon further straining, but their growth rapidly 
becomes arrested Q- Instead of forming new voids, ad- 
ditional material is drawn out of the uncavitated poly- 
mer, and stable craze growth occurs. In our simulations, 
growth continues until all material in the simulation box 
is converted into the craze. 



A. Images of crazes 

A good impression of the crazing process can be ob- 
tained by inspecting the snapshots of the simulation cell 
shown in Figs. 0] - HO Each slice has a lateral width of 
64 a, and three different strains are shown. In all im- 
ages, the chain length N — 512. Previous studies [9] had 
shown that N has to be twice the entanglement length or 
greater in order to form stable crazes. For shorter chains, 
the material cavitates, but then rapidly fails due to chain 
pullout. In the following, we only consider chains with 
N > 2N e . 

Note first that in all cases, there is a sharp inter- 
face between dense polymer and crazed material. This 
narrow 'active zone' is one of the key features of craze 
phenomenology found in experiment. In the craze, the 
polymer chains have merged into fibrils that are strongly 
aligned. However, the structure is quite complicated, as 
there are many lateral connections between fibers. 

One can also observe that the fine structure of the 
crazes in the three sequences varies greatly. Fig. 0] with 
flexible chains at the low temperature of T = O.luo/ks 
and the weak adhesive interaction (cutoff distance r c = 
1.5a) shows many thin fibrils, whereas the fibrils in 
Fig. [fj] at the higher temperature of T = 0.3uo/^s and 
the stronger adhesive interaction r c = 2.2a are much 
thicker in diameter. These trends are not surprising, be- 
cause increased chain mobility at higher temperatures 
and stronger adhesive interactions should drive the sys- 
tem to larger fibril diameters, which minimize the surface 
area. 



B. The drawing process and stress-strain curves 

A second characteristic feature of craze growth is that 
deformation occurs at a constant plateau or drawing 
stress S. This plateau can be easily identified in the 
stress-strain curves shown in Fig. [7] The curves can be 
separated into three different regimes. In regime I, the 
stress rises to a peak of ~ 2.6uo/a 3 and then drops when 
the polymer yields by cavitation. Following cavitation, 
the stress rapidly relaxes and remains at the plateau 
value S in regime II, the growth regime. Regime II is 
much shorter in the semiflexible case Fig. Efb) than in 
the flexible case Fig. [7fa) (note different lateral scales). 



FIG. 4: Three snapshots of craze growth for flexible chains 
with T = O.lwo/fcs and r c = 1.5a. The total system contains 
262144 beads, but only slices of thickness 10a normal to the 
page are shown in order to resolve the fine structure. The 
lateral dimension of each slice is 64a and the vertical direction 
is to scale. Each dot represents one Lennard- Jones bead. 




FIG. 5: Three snapshots of craze growth for semiflexible 
chains with T = 0.1-ao/fcs, r c = 1.5a, and 262144 beads. 



Regime II ends when the strain L z /L reaches the exten- 
sion ratio A. At this point, all the material in the simula- 
tion cell has been converted into the craze, and additional 
deformation strains the entire craze uniformly. As a con- 
sequence, the stress rises again in regime III. This regime 
finally ends in catastrophic failure either through chain 
disentanglement or chain scission (see Section IVlTl) . 

Note first that neither the peak stress at cavitation 
nor the value of S depends on the chain length N . The 
curves for different N in Fig. [7] only split apart after 



FIG. 6: Three snapshots of craze growth for semiflexible 
chains with T = 0.3uo/fcs, r c = 2.2a, and 262144 beads. 



completion of craze growth when L z /L reaches A and the 
entire craze is strained. Baljon and Robbins [9j showed 
that the peak stress remained constant for much shorter 
chains, but that regime II only appeared when N was 
2N e or longer. Another important fact to note is that 
S is independent of system size. For example, values of 
S in systems ranging between 32768 and 1048576 beads 
are the same within a few percent. The biggest change 
with increasing system size is that temporal fluctuations 
in S decrease. 

In Fig.|Hta), we analyze trends of S with T and r c . The 
drawing stress decreases linearly with increasing temper- 
ature and increases with increasing adhesive interactions 
(i.e. increasing r c ). Fig. (b) shows that S varies log- 
arithmically with the widening velocity v over two or- 
ders of magnitude, which is indicative of a thermally ac- 
tivated process. For the subsequent figures, we choose 
v = 0.06a/rLj, which is at the upper end of the loga- 
rithmic regime || . Similar behavior is also found for the 
shear yield stress of glassy polymers [13, ■ 

C. Crazing under plane stress conditions 

The results of Section lTVI show that cavitation only oc- 
curs when all three principal stresses are tensile. Many 
experimental crazes grow in a thin film geometry under 
plane stress conditions. However, in these experiments 
the craze is often prenucleated or nucleates near a defect 
|12|. This situation can also be mimicked in our simu- 
lations. To this end, the periodic boundary conditions 
in the x direction were replaced with free boundaries, so 
that the solid is free to relax in that direction. Initial 
failure is now nucleated by placing 1000 purely repulsive 
LJ beads in the center plane of the simulation cell located 
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FIG. 7: Stress cr zz in the widening direction during craze 
growth at T = O.ltto/fcs, r c = 1.5a for (a) flexible and (b) 
semiflexible chains of length TV = 128, AT = 256, TV = 384, 
and TV = 512. Three characteristic regimes of (I) cavity nu- 
cleation, (II) craze growth and (III) craze failure are also in- 
dicated. The two perpendicular stress components a xx and 
a yy also peak at cavitation (see text), but then rapidly drop 
to zero. Qualitatively identical curves are obtained at other 
values of T and r c . 



at z = L z /2 [4j|. This weakens the solid locally and con- 
strains the location of initial failure, while not affecting 
subsequent craze growth. 

Fig. shows three snapshots of a craze in this geome- 
try. As in experiments, necking is observed at the craze- 
bulk interface. Although er xx vanishes in the rest of the 
film, the neck produces strong tensile stresses in all three 
directions in the active zone. The craze grows in the same 
fashion as in the simulations with 3D periodic boundary 
conditions. Since the latter yield better statistics for the 
craze structure, we have focused on this methodology for 
our analysis. 



D. The extension ratio 

The extension ratio A can be calculated from the av- 
erage densities of crazed and uncrazed material. Fig. I1UI 
shows how the density drops from the initial value pi to 
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FIG. 8: (a) Trends of S with T and r c at v = 0.06a/r LJ 
for flexible (■) and semiflexible (a) chains and r c = 1.5a 
(lower curves) and r c = 2.2a (upper curves), (b) Velocity 
dependence of S for flexible chains at T = O.luo/kB. The 
straight line is a fit to a logarithmic velocity dependence, S = 
1.085 uo/kB +0.048 uo/kB lnv. Uncertainties are comparable 
to symbol sizes. 



Pf in the craze. As can be seen, pf is higher for the semi- 
flexible chains, which have a smaller value of TV e w 30. 
Remarkably, we find that A is a function of N e only and 
decreases with decreasing N e . For instance, while in- 
creasing T and r c produces dramatic coarsening of the 
fibril structure in Fig. Irrelative to Fig.0J A is unchanged. 
We obtain values of Aa = 6.0 ± 0.6 and A s fl = 3.5 ± 0.3 
independent of TV , T, and adhesive interaction strength. 

In order to understand the dependence of the macro- 
scopic quantity A on TV e , we analyze the structural 
changes in the polymer glass during deformation on a mi- 
croscopic level (see also ref. Figure UTT a) shows the 
average final position of beads in the completely evolved 
craze as a function of their initial positions along the di- 
rection of the expansion (z-axis). The average was taken 
over all beads with initial heights in a bin of width la. 
Although the strain rate is strongly localized during the 
craze process, the ultimate displacement profile is linear, 
z f = Xzi. 

To measure deviations from a purely affine (uniform) 
deformation, we evaluated the rms variation Sz in Zf for 
beads in each bin. This quantity is indicated by error 
bars in Fig. Illf a). Note that the variation in each bin is 
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FIG. 9: Cross-sections through a craze with a free interface 
at T = O.luo/fcs, r c = 1.5a, and 262144 beads. Periodic 
boundary conditions were maintained in the direction into 
the plane. The location of initial cavitation was constrained 
by placing repulsive beads in the center plane at z = L z /2. 
The lateral dimension is 47 a and the vertical dimension is to 
scale. 




20 40 60 80 100 120 

z/a 

FIG. 10: Density profile through the active zone for crazes 
with flexible chains (N e « 70) and semiflexible chains (N e ~ 
30). Horizontal lines indicate the average density in the craze 
for the two cases. 



very reproducible. We find that Sz is nearly independent 
of T and r c and has values on the order of 19a and 9a for 
flexible and semiflexible chains, respectively. 

Since no strain is applied in the perpendicular x and y 
directions, one would assume that there is on average no 
displacement in these directions. That this is indeed the 
case is shown in Fig. ^] which repeats the analysis of 
Fig. ^] for the x-direction. Average final bead positions 
are identical to initial positions, but there are lateral vari- 




FIG. 11: (a) Final bead heights Zf as a function of initial 
heights Zi for flexible (large slope) and semiflexible (small 
slope) chains (T — O.le/fcs, r c — 1.5a). Averages were cal- 
culated over z-intervals of width a. Straight lines have slope 
A = 5.9 and A = 3.5, respectively. Error bars represent a 
standard deviation from the averages in each layer and are on 
the order of 19a (flexible) and 9a (semiflexible). (b) Square 
of the height change Az as a function of the number of cova- 
lent bonds AiV between a bead and the chain center. Dashed 
straight lines have slope \ 2 l p lo/3 with A from (a). Deviations 
from the RW scaling occur in the vicinity of the chain ends 
(not shown). Other systems at different T, r c and N show the 
same results. 



ations Sx that are indicated by error bars. These lateral 
displacements allow chains to gather in fibrils at the ini- 
tial density to minimize surface area. Unlike the verti- 
cal displacements Sz, these lateral displacements depend 
strongly on T and r c . For example, Sx ~ 2.5a for the fine 
structure shown in Fig.0J where many thin fibrils can be 
seen, while Sx ~ 5.6a for the much coarser structure of 
Fig. El In general, Sx correlates with the spacing between 
fibrils as discussed in Sect ion IVTl and is less than dj. Kru- 
penkin and Fredrickson suggested that di provides an 
upper bound for the lateral chain deformations. 

We now examine changes in the conformation of in- 
dividual chains. In the initial state, the polymer chains 
exhibit an ideal random walk (RW) structure inherited 
from the melt. The average end-to-end vector (R 2 ) thus 
scales with the number of covalent bonds connecting two 
beads AN as (R 2 ) = I p IqAN. The component along 
each direction is 1/3 of that value since the initial state 
is isotropic. Fig. ll2f b) shows this initial scaling behavior 
for (Ax 2 ) (dashed line) and that (Ax 2 ) is not affected by 
crazing (solid line). 

After an affine deformation by A along z, one would 
have an anisotropic RW with no change in Ax or Ay, 
but (Az 2 ) = X 2 l p l AN/3. Fig. ITlTbl shows the actual 
behavior (solid lines) of (Az 2 ) in the craze. At large 
scales, it exhibits the expected scaling for an affine de- 
formation (dashed lines). However, the separation be- 
tween beads is fixed by the length of the covalent bonds, 
so the deformation of individual polymers along z can- 
not be purely affine. At small scales, the linear scaling 
behavior of (Az 2 ) crosses over into a quadratic behav- 
ior, which indicates that the polymer has been pulled 
taut on this scale. The typical number of beads in 
such a straight segment N st can be calculated by let- 
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Xi/a AN 

FIG. 12: Analysis of bead positions analogous to the previ- 
ous figure (same systems), but for the x-positions. No strain 
is applied in this direction, and the straight lines in panel (a) 
have slope one. The curves for the semiflexible chains in (a) 
were displaced vertically upward by 10a to avoid overlap. Er- 
ror bars represent a standard deviation Sx from the averages, 
and are on the order of 3.7a (flexible) and 2.2a (semiflexible). 
(b) Bead displacements as a function of distance from the 
center in bond lengths, AiV, along the chain. Dashed lines 
have slope l p lo/S. 



ting (7V st /o) 2 = (Az 2 ) = \ 2 lpl N st /3 at the crossover 
point, which yields N st = X 2 I p /31q. Inserting the ob- 
served values of A, l p , and Iq, we arrive at N^ t = 21 ± 4 
and = 12 ± 2, respectively. These values are com- 
parable to the values of Sz found in Fig. Illf a) . On this 
scale, the deformation is non-affine. 

The length of taut sections can also be determined by 
direct analysis of the chain geometry. To this end, we 
calculate the angle between every covalent bond and the 
z-axis and label a bond as pointing up (down) if the an- 
gle is within 45° of the z (-z) axis. We then count the 
number N st of consecutive up (down) steps. The prob- 
ability P(N st ) of finding a straight segment containing 
N st steps is shown in Fig. I13f a). For both flexible and 
semiflexible chains, the distribution develops an expo- 
nential tail. Like A, this tail is independent of N, T, and 
r c . The characteristic length scales that arise from these 
tails are N^t ~ 21 and iVjf ~ 13, in good agreement with 
the prediction from the RW argument. Fig. I13f b) shows 
that very similar length scales arise from an equivalent 
analysis of the decay of the correlation function for the 
z-component of successive bonds. 

In section III B II we introduced the standard scaling 
argument Eq. J3J) that relates extension ratio and entan- 
glement length, which has been verified experimentally 
with great success. In our cases, it predicts A„ lax = 6.5 
and A^ ax = 3.5, which agree with the observed values of 
A. However, the argument was motivated by the idea that 
segments between entanglements become fully stretched 
and thus appears to be at odds with the finding of an av- 
erage straight segment length of only N e /3 rather than 
N e . This discrepancy is resolved by realizing that since 
the deformation is uniaxial, only the projection of di onto 
the z-axis, diCos(Q z ), is expanded, where Z is the an- 
gle between di and the z-axis. The average projection 
is thus only l/\/3 of the total length. Indeed, it was 




FIG. 13: (a) Probability distribution of straight segments of 
length N s t for flexible and semiflexible chains. Thick lines cor- 
respond to simulations at T — O.le/fcs, r c = 1.5a, N = 512 
with 1048576 beads. Dotted lines were obtained at T = 
0.3e/fc s , r c = 2.2a, N = 512 with 262144 beads, and long 
dashed lines correspond to T = O.le/fcs, r c = 1.5a, N = 
256 with 262144 beads. The straight lines show fits to 
exp(— Nst/N^t ). (b) z-component of the bond-bond corre- 
lation function for the same systems. Thin solid lines show 
exponential fits with the indicated decay lengths. 



noted already in earlier work [1J| that, due to this ge- 
ometric factor, A should be y/3X max for fully stretched 
chains. However, this result is little cited since A ~ A max 
in many systems and, until our work, there was no reason 
to expect the length of straight segments to be N e /3. 

The emergence of the length scale N e /3 is a conse- 
quence of the random nature of the entanglement mesh. 
Clearly, all strands would be expanded simultaneously by 
the same factor in a regular mesh as reported in a simu- 
lation study by Stevens [27]] ■ By contrast, in the polymer 
glass only the segments that are initially aligned with the 
stretching direction become fully stretched. These fully 
stretched segments are able to prevent further extension, 
because the entanglements act like chemical crosslinks. 
Barsky and Robbins have confirmed the equivalence be- 
tween entanglements and crosslinks by adding permanent 
crosslinks randomly to the system (37|. The length be- 
tween constraints then decreases from N e , and A max de- 
creases accordingly. They found A ~ A max in all cases 
and that the average stretched length N st remains at 1 /3 
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TABLE I: Dissipation during craze growth and covalent con- 
tribution to the crazing stress S for several different systems 
of size 262144. 

T r c N SQ/SW % cov stress 



fl. 


0.1 


1.5 


256 


0.88 


87 


fl. 


0.1 


1.5 


512 


0.88 


88 


sfl. 


0.1 


1.5 


256 


0.71 


95 


sfi. 


0.1 


1.5 


512 


0.67 


97 


fl. 


0.1 


2.2 


512 


0.92 


69 


sfi. 


0.1 


2.2 


512 


0.71 


75 


fl. 


0.3 


2.2 


512 


0.87 


61 


sfi. 


0.3 


2.2 


512 


0.78 


67 



of the distance between constraints. 

The success of the scaling argument Eq. (J2J and the 
constancy of the extension ratio imply that there is no ap- 
preciable loss of entanglements in our simulations during 
craze growth. Chains do not disentangle once N > 2N e , 
and chain scission (see also section lYlIfl is not observed 
during growth for any choice of parameters in our model. 



E. Energy dissipation and stress transfer during 
crazing 

The work done in transforming a volume dV of poly- 
mer into a craze is SW = S(X — \)dV . This work can 
either increase the potential energy dU or be dissipated 
as heat SQ. The division between energy and heat is hard 
to determine experimentally, but simulations with short 
chains found that both contributions were substantial @ . 
We have measured SW and the energy change directly in 
our simulations and calculated SQ using the first law of 
thermodynamics: SQ = dU — SW. dU can be calculated 
directly from the bead positions and interaction poten- 
tials. Table [I] shows the fraction SQ/SW of dissipated 
total work for a number of large systems. In all 
large percentage, ~ 80%, of the total work is dissipated, 
and only ~ 20% is stored as potential energy. Since the 
craze drawing stress varies logarithmically with velocity 
(see Fig. ISfb)), these percentages could change with ve- 
locity. However, we find that dU also decreases with 
decreasing velocity, and there is no measurable change 
in the percentage of work converted to heat over at least 
two orders of magnitude in velocity. 

Stress in the craze can also be partitioned into two 
components that originate either from van der Waals 
(LJ) interactions Eq. © or from covalent interactions 
Eq. Q. The two contributions are very different in 
the uncrazed and crazed material. In the undeformed 
polymer, the tensile stress is mainly carried by the vdW 
bonds. This stress is transferred in large part to the 
covalent bonds as they pass through the active zone. 
Evidence for this is provided in Fig. 1141 which shows 
the covalent and LJ contributions to the total stress as 
a function of height in the widening direction. Panel 
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FIG. 14: (a) Density profile through a craze simulation at 
T = O.luo/kB, r c = 1.5a with flexible chains, (b) vdW stress 
(solid line), covalent stress (dashed line) and total stress S 
(thick line) as a function of position along z. The kinetic 
contribution to the stress is split evenly between the covalent 
and vdW stress here and in Table Q 



(a) displays the density profile in order to identify dense 
polymer regions (high density) and craze regions (low 
density). In Panel (b) one observes that in the dense re- 
gion all the tension is carried by the van der Waals bonds 
and the covalent bonds are under slight compression. In 
the craze, between 60 — 95% of the total stress (see Table 
HJ) is carried by the covalent bonds, and the van der Waals 
bonds only contribute a small fraction. 



F. Problems with the capillary models 

The results presented so far reveal serious difficulties 
with the surface tension models discussed in the Intro- 
duction. The first evidence of this comes from the obser- 
vation that S is independent of system size. In our small- 
est simulations of lateral width 32a the simulation box 
only contains a few fibrils at T = OAu^/kg, r c — 2.2a. If 
S were controlled by Eq. Q, one should expect that the 
simulation box would need to contain a statistically sig- 
nificant number of fibrils of spacing (-Do) f° r S to reach 
its steady state value. However, the value of S does not 
fluctuate as the lateral dimensions are increased to 64a 
or 128a. 

The second and more severe problem concerns the dis- 
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tribution of stress in the craze. The surface tension model 
assumes that all the stress is carried at the interfaces of 
fibrils or in viscous stress in the active zone. However, 
Tableland Fig. 1 141 show that almost all the stress in the 
fibrils is carried by the covalent bonds, while the surface 
tension is entirely associated with broken vdW bonds and 
small entropic contributions. 

In the following, we make an alternative proposition 
to relate craze microstructure and drawing stress. This 
proposition is based on the observation that the values of 
S and A obtained from our simulations obey the equality 

SflAfl = SsflAsfl = S (T, r c , v). (12) 

This can be verified for each T and r c using A a = 6.0, 
A s fl =3.6 and values of S from Fig. EI Since the fraction 
of area occupied by the fibrils is 1/A, So is the local stress 
within the fibrils. 

It is perhaps surprising that the value of local stress 
needed to draw fibrils is independent of the entangle- 
ment length. In order to further test Eq. 1|12[> . it would 
be desirable to consider additional values of N e and thus 
A. Unfortunately, reliable values for N e exist for only a 
few values of b. It is, however, not necessary to know 
N e for the present purpose, since both S and A can 
be measured directly from the craze simulation. More- 
over the entanglement length should only depend upon 
the chain statistics 36], and glassy states with arbitrary 
statistics can easily be created. We confirmed that sim- 
ulations with the same persistence length in the unde- 
formed glass gave the same values of A and S indepen- 
dent of whether the bond-bending potential (Eq. (|SJ)) was 
included. Fig. llSf al compares stress strain curves for 
b = Ouq at four values of the initial persistence length. In- 
creasing l p lowers A and increases S. However, Fig JlSf b) 
shows that all curves can be collapsed if a zz is scaled by 
A and the extension by A — 1. This confirms that So is 
the stress that controls craze growth and only depends on 
the van der Waals interactions. An experimental version 
of this test would be difficult, since it is hard to change 
N e for real polymers without changing the chemistry as 
well. 

It is interesting to compare values of So to the stresses 
required for shear yielding and cavitation, which arc also 
independent of chain statistics. Table ITT1 compares these 
three stresses for two ranges of the LJ potential. Here 
b = Quo, but the bond-bending potential has little effect 
on the values. The three stresses are clearly correlated, 
decreasing with increasing temperature and decreasing 
r c . The local fibril stress is always about twice the cav- 
itation stress and ranges from 7 to 12 times To. The 
implication is that local stress for drawing material into 
fibrils So is related to the bulk yield stresses, but it is 
difficult to determine the relative role of shear and cavity 
growth. It is interesting to note that the experimental 
values of S /r for PS and PMMA are about 5 0IH, 
and it would be useful to have values of this ratio for 
other polymers. A comparison to p cav would also be in- 
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FIG. 15: (a) Stress-strain curves for crazes with l p — 1.65a, 
2.2a, 2.7a, and 3.3a in order of increasing height at T = 
O.lwo/fcs, t c = l-5a and b = Otto- The corresponding val- 
ues of A are given in Table |TTT| (b) Rescaling of the same 
data in the form a zz \ versus (L z /L — 1)/(A— 1). All curves 
collapse onto a common plateau So = SX ~ 4.5uo/fce of the 
same length. 



TABLE II: Values of the shear yield stress to, the yield stress 
for cavitation p cav and the local fibril stress So as a function 
of temperature. Stresses are in units of Mo /a 3 . Uncertainties 
are ±0.02 in to, and about 10% in the other quantities. 



Tk B /u 


r c — 1.5a 
0.3 0.1 0.01 


r c = 2.2a 
0.3 0.1 0.01 


To 

Pcav 

So 


0.23 0.49 0.72 
1.2 2.7 3.0 
2.9 4.5 5.2 


0.45 0.64 0.83 
3.0 4.8 5.0 
5.9 8.0 



teresting, but its value is sensitive to system size, strain 
rate and inhomogeneities and it is difficult to measure 
experimentally. 

Given that previous results for N e [2^ |2|| |2j| are con- 
sistent with values inferred from the extension ratio (Eq. 
|2J), our results for A as a function of l p allow a rapid 
estimation of N e . Table II I II presents results for a wide 
range of l p and shows that the product N e l p is constant 
within our errorbars. Fetters et al. 36] have presented 
a model for the relation between chain statistics and N e 
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TABLE III: Measured values of A, with uncertainties, as a 
function of l p and the corresponding range of values of N e and 
l p N e inferred from Eq. [5] Runs were made at T = O.luo/fcs, 
r c — 1.5a and b = Ouq. 



A 



1.65 6.0 ± 0.6 50-76 83-124 

2.2 4.5 ±0.5 37-57 81-127 
2.7 3.5 ±0.3 29-41 78-111 

3.3 3.0 ± 0.3 25-37 83-122 
5.55 2.0 ±0.2 19-28 105-155 



that predicts 



N K 



N 3 

WY 



(13) 



where (R 2 ) = l p l$N denotes the average end-to-end vec- 
tor of the polymer chain and the proportionality constant 
only depends upon density. This implies N e oc l~ 3 , while 
our data is clearly consistent with a simple inverse re- 
lation N e on lp 1 . Equation l|13|) describes many experi- 
mental polymers, but it is difficult to change l p without 
changing all the other parameters in the equation. The 
flexible model (b = Ouq) is known to be quantitatively 
inconsistent with Equation 113|) 25] , which has been one 
motivation for studies of more rigid models. It would 
be interesting to have additional values of l p from melt 
simulations to test whether the inverse relation between 
N e and l p found here holds more generally, and, if so, to 
understand its origin. 



G. Width of the active zone 

At the interface between dense polymer and craze, 
polymer chains are locally mobilized and brought into 
the new fibril structure. The region in which this motion 
takes place is called the 'active zone'. In Fig.|2the height 
of the active zone h was defined as the distance between 
undeformed polymer layer and the void ceiling, and this 
layer was assumed to behave like a strain-softened fluid. 
The main drop in density should occur over the height of 
the void ceiling of order ~ (Do)/ 2 ^ n this section, 

we compare this simple picture to our simulations. 

Fig.llOlshows typical results for the density profile near 
the craze boundary. For both flexible and semiflexible 
polymers, the density drops over a region of width ~ 20a. 
The average strain rate e must be localized in the same 
region, since e = —d\np/dt. If the active zone advances 
at velocity v then e = ±vdhip/dz (the sign depends 
on whether the top or bottom interface is growing). In 
Fig. Iltjf a) we present e as a function of position along 
z. Averages were taken over layers of height la. The 
curves shown correspond to the 4 different values of l p 
used in the creation of the polymer chains. As discussed 
above, this varies the entanglement length N e . Curves 
were shifted by z$, which corresponds to the center of 
the peak. At this point, the density is roughly halfway 
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FIG. 16: (a) Strain rate e as a function of distance from the 
location of the onset of density drop zo for 4 values of N e 
(see text) . Positive values of z — zo measure the distance into 
the crazed region. The length of the dashes increases with 
decreasing N e . (b) Standard deviation a(Az) from the aver- 
age displacements along z in a given layer over a time interval 
of 75tlj. ct(Az) decays exponentially into the craze and the 
straight lines show fits to this decay with characteristic length 
scales 48a, 35a, 27a, and 24a. 



between pi and pf. Curves for very different N e essen- 
tially overlap with a width of ~ 10a at half maximum 
and ~ 20a at the base. The latter corresponds to the 
range of rapid density change in Fig. QJJ 

The width of the region over which beads are mo- 
bilized was determined from the relative diffusion as a 
function of height. Fig. Iltjfl j) shows the standard devia- 
tion of displacements in the lengthening direction cr(Az). 
The curves peak at the same location as the curves in 
Fig. Iltjf al. but are more asymmetric. In the direction of 
the dense polymer glass ct(Az) and epsilon fall to zero 
over a comparable range. By contrast, a(Az) shows a 
long exponential tail into the craze with a characteristic 
decay length Iaz that varies with entanglement length. 
The fit values of Iaz indicate that there is a definite trend 
to larger values as N e increases, and Iaz tends to be 
somewhat smaller than N e . This result is not surprising, 
because N e is the longest length scale over which parti- 
cle mobilization should occur. Standard deviations of the 
lateral displacements Aa; and Ay are smaller, but show 
essentially the same decay lengths. 
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The above analysis indicates that while the mobility of 
the beads is constrained by entanglements, the regions of 
localization of strain rate and the density drop are related 
to the craze microstructure (see Section lVl*|) . Typical val- 
ues for (D) and (-Do) are given in Table llVl From this, 
(Do) /2 ~ 10a, which compares well to the width of the 
strain localization peak at half maximum. Experimen- 
tally, the width of the active zone has been measured by 
a gold decoration technique 0] . It was concluded that it 
lies between (D) and (Do), which agrees with our results. 



VI. MICROSTRUCTURE OF CRAZES 

Another fascinating aspect of crazes is their compli- 
cated microstructure. Figs. 01 - El give an impression of 
the range of length scales appearing in the voided fib- 
ril network. Clearly, the picture of cylindrical fibers and 
void fingers (Fig. [3J is an oversimplification. It is nev- 
ertheless helpful to build more realistic models starting 
from this simple scenario. 



A. Structure factor 

Experimentally, the standard measurement of the 
craze microstructure is done via scattering experiments. 
The scattering intensity in these measurements is pro- 
portional to the structure factor 



5(k) 




ik(rj-r fc ) 



(14) 



times the form factor for the monomers. Since the craze 
structure has azimuthal symmetry, one decomposes the 
wavevector k into components parallel and perpendic- 
ular to the fibrils. Contour plots of S(k±,k^) for two 
crazes are shown in Fig ll7l The microstructure was var- 
ied by changing the cutoff length r c and chain flexibility. 
Both patterns are asymmetric with the intensity decay- 
ing much faster in the direction parallel to the fibrils than 
perpendicular to them. 

Most experimental setups integrate over fcy using slit 
collimation and measure the integrated structure factor 
S(k 1 _) = f^dk^Sik). In Fig. HHJa), we plot S(k x ) as 
a function of the magnitude of k± for the same systems 
shown in Fig. El At large wavevectors, the curves rise 
to a peak at ~ 2ir/a (not shown), which corresponds to 
the typical separation of two beads. This length scale is 
so short that it is usually not resolved in typical experi- 
mental scattering plots shown in e. g. [2lj . The charac- 
teristic feature of S(k±) is found at smaller k± in form 
of a power-law regime with exponent —3. The extent of 
this scaling regime is bound at large wavevectors by the 
small scale cutoff provided by the interparticle spacing 
and at small wavevectors by the distance to the next fib- 
ril. The power-law regime is more pronounced for the 
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FIG. 17: Contour plots of the structure factor of crazes with 
(a) flexible chains, r c = 1.5a and (b) semiflexible chains, 
r c — 2.2a. The temperature was T = O.luo/kB and the 
systems contained 1048576 beads. Colors range from black 
(high intensity) to white (low intensity). 



crazes shown in Fig. I18f b) , which were created at higher 
temperatures where (Dq) is bigger. 



B. Interpretation of scattering data 

The traditional interpretation of craze scattering data 
begins with idealizing a craze fibril as a straight cylinder 
of diameter D and length I along z. The scattering in- 
tensity is then proportional to the squared magnitude of 
the form factor for such a cylinder 



jquar 



F(k x ) 



ApeiirD 2 ! 1 / 2 JxiirDkx) 



(15) 



where Ap e i denotes the electron density and J\ is the 
first-order Bessel function. Due to the asymptotic be- 
havior of Ji(x) = (2/7rx)V 2 cos[a: - 3tt/4] + O(x^) 
for large arguments, the scattering intensity of a single 
cylinder will exhibit an oscillating power-law behavior, 
\F(k±)\ 2 oc fc7_ 3 - This is also called Porod scattering. 

In general the fibrils do not have a single diameter, but 
rather a diameter distribution P(D). One can introduce 
an average scattering intensity 



/o(*l) = (F(k ± ) 2 ) 



dDP{D)F{k A 



(16) 



by averaging the form factor over the diameter distri- 
bution and neglecting correlations between fibrils. The 
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FIG. 18: Integrated structure factor of crazes (a) at T = 
O.lMo/fes, flexible chains, r c — 1.5a (solid) and semiflexible 
r c — 2.2a (dashed) and (b) at T = 0.3ito/A;.B, flexible (solid) 
and semiflexible (dashed) chains with r c = 2.2a. The straight 
lines have slope -3. System sizes were 1048576 beads in (a) 
and 262144 beads in (b). 



main effect is to smooth out the oscillations so that a 
straight power-law tail results. The average diameter 
(D) = J DP(D)dD can be obtained via a Porod analysis, 
in which one determines the prefactor a to the power-law 
tail, S(k±) = afcj 3 . This can be related to (D) [3113 
through, 



(D) = 



tt 3 (1 - \/X)a 



(17) 



where Q = J dk±2irk±S{k±) is a scattering invariant. 
Values for (D) obtained from this formula are collected 
in Table ITVl 

The craze fibrils do not all have the same distance from 
each other, but have in general varying distances that can 
be described by a radial distribution function g(r). This 
will lead to interference effects in the scattering intensity 
that can be described by an interference function j(k±) = 
I(k±)/In(k±) — 1. The interference function is related to 
g{r) by [2l| 



l(r) 



1 



ir(D} 2 \ 



2irk ± j(k ± )J (2irrk ± )dk ± . (18) 



4'* ' 



As a result, the power-law tail will be modified at small 



FIG. 19: Typical crosssection through a craze from flexible 
chains, r c — 1.5a, T = O.lwo/fcs ■ The lateral dimension 
is 128 a. Fibrils appears as clusters of varying size, and the 
distributions in Fie. l20l are calculated form these crosssections 
(see text). 



wavevectors. In particular, the first peak in g(r), corre- 
sponding to a typical fibril separation {Da), should trans- 
late into a maximum in I{k±). 

The Porod scattering law |F(fcj_)| 2 cx kj 3 is well con- 
firmed at higher temperatures in Fig. I18f b). while the 
power-law regime is shorter at lower temperatures. Vi- 
sual inspection of the craze images suggests that the 
straight cylinder approximation is not so well justified 
in this case. Fibrils branch more often and intersect the 
z-axis at varying angles. At the higher temperature, the 
chains are more mobile and can align more easily, but 
they are still not ideal cylinders. We note furthermore 
that the curves shown in the log-log plot of Fig. ^3] do 
not exhibit a clear maximum (a maximum would be more 
easily identifiable in a linear plot normally used for ex- 
perimental results). This would suggest that the ordering 
of fibrils is mostly random without a clear characteristic 
separation. However, our statistics are limited by the 
system size at these large length scales. 



C. Distributions of fibril diameter and spacing 
from real space analysis 

In a previous analysis of scattering data [2(| for 
polystyrene and polycarbonate crazes, Eqs. JTSJ-lJTSJl 
were used to extract the diameter distribution P{D) and 
g{r) by means of a detailed fitting procedure. The craze 
images shown above suggest that P(D) is rather broad, 
and wide distributions for P{D) were found for both ma- 
terials, with a significant increase in breadth for polycar- 
bonate [2l| . 

Here, we access these distributions by direct geomet- 
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FIG. 20: Distribution of fibril diameters from analysis of con- 
nected clusters (see text) for the same crazes as in Fig. I18f a~l 
and (b). Solid lines refer to flexible chains, dashed lines to 
semiflexible chains. 



rical analysis of the bead positions. To this end, we bin 
the particle positions onto a square grid with gridsize la 
normal to the widening direction and take lateral cross- 
sections of height 1.5a. As illustrated in Fig. El a fibril 
now appears as a 2-dimensional connected cluster, whose 
area A is taken to be the s um of t he areas of the occupied 
squares. We define D = y/4A/n. 

Fig. 1201 shows the resulting distributions of D. As ex- 
pected, P(D) is very broad. The distributions for flexible 
and semiflexible chains are very similar at small D for a 
fixed temperature, but differ for larger D. The tail of the 
distributions could be fitted to an exponential function, 
but our statistics are too small for a conclusive state- 
ment. Mean values of the diameter are given in Table 
dVI together with the standard deviations of the distri- 
butions. The large values of the latter suggest that (D) 
has to be used with care when describing the craze mi- 
crostructure. Previously, Baljon and Robbins reported 
similar values of (D) — 7a, cr(D) = 11a for flexible chains 
at T = 0.3uo/fcs, r c = 1.5a @|. 

The value of (D) obtained from the Porod analysis is 
always larger than the value from the cluster (real-space) 
analysis. Both values rise with increasing adhesive inter- 
action and increasing temperature as expected. Taking 
a ~ 0.8nm, diameters in our model crazes would corre- 



TABLE IV: Structural parameters of model crazes. Size refers 
to the total number of beads in the simulation. For the fibril 
spacing, results for (D) from both the scattering analysis and 
the cluster analysis (see text) are shown. The rms variation 
a(D)/a was obtained from cluster analysis. 

T r c size {D)/a (D)/a a(D)/a (D )/a 
scatt. cluster 



sfl. 


0.3 


1.5 


128 s 


6.1 


5.1 


9.4 


18.4 


sfl. 


0.3 


2.2 


128 s 


11.3 


6.8 


13.1 


25.0 


fl. 


0.1 


1.5 


64 x 128 2 


4.7 


4.2 


4.1 


14.3 


sfl. 


0.1 


1.5 


64 x 128 2 


4.8 


4.0 


5.5 


12.8 


fl. 


0.1 


2.2 


64 x 128 2 


8.4 


5.5 


7.0 


22.3 


sfl. 


0.1 


2.2 


64 x 128 2 


8.2 


5.3 


9.5 


19.8 


fl. 


0.3 


2.2 


64 3 


12.6 


7.9 


8.7 


30.7 


sfl. 


0.3 


2.2 


64 3 


11.1 


6.5 


10.4 


23.5 
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FIG. 21: Radial distribution function g(r) from analysis of 
connected clusters (see text). The systems shown are (a): T = 
O.lwo/fes, r c = 1.5a, flexible (♦) and r c = 2.2a, semiflexible 
(*), (b): T = 0.3ii /fc s , r c = 1.5a (♦) and r c = 2.2a (*) (both 
semiflexible) . 



spond to a range of (D) = 1.4 — 6.3nm, which is at the 
small end of the experimental range. The reason is that 
an artificially small value of the surface tension and a 
high widening velocity v are used here. Both lower (D), 
which allows us to use smaller system sizes. 

An estimate for the mean fibril spacing (Do) can be 
obtained by equating the area per fibril, ttDq/A, to the 
inverse areal density 1/n, i.e. (Do) = 2^/1/mr. The areal 
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density was obtained by counting the number of separate 
fibrils per cross-section. Values for (D ) are also given in 
Table lfVl and translate into a range between 10 nm and 25 
nm. The higher numbers are comparable to experiment 
and are obtained with r c = 2.2a, which produces more 
realistic surface energies. 

In order to obtain the radial distribution function of 
the fibrils, we continued the analysis described above and 
calculated the center of mass for each 2D cluster. The 
positions given by this procedure were used to calculate 
g(r) in Fig. 1211 In general, these functions have very 
little structure. There is a size exclusion minimum at 
the origin, and the curves have a weak first maximum 
around 10a. As the fibrils become thicker, the location 
of the maximum shifts to larger values. Qualitatively 
similar curves were obtained from experiment |2l| , which 
confirms the basically random nature of fibril positions. 
The height of the maximum is too small to be reflected 
in the scattering intensity. 




5 10 15 20 25 

N/N e . 

FIG. 22: Saturation of fibril breaking stresses in systems of 
size 262144 beads, T = O.luo/fcs. 5*00 denotes the maximum 
saturation stress in the limit of very long chains. The solid line 
is 1 — cxp (—N/4N e + 1/2). Squares indicate flexible chains 
(N e ~ 64) and triangles indicate semifiexible chains (N e ~ 
32). 



VII. STRESS DISTRIBUTION AND CRAZE 
BREAKDOWN 

In regime III of the stress-strain curve of Fig. [7| the 
entire volume of the simulation cell has been converted 
to a craze. Elongation past the extension ratio causes 
uniform straining of the craze and eventually leads to 
craze failure. Studies of this regime are directly relevant 
to crack propagation in glassy polymers (Fig. [y. The 
stress in the craze region rises from S at the active zone 
to a maximum value S max at the crack tip. The elastic 
properties of the craze determine the rate at which the 
stress rises with distance, and S'max determines how large 
the craze region can become before the crack propagates 
[||. These properties were recently obtained from MD 
simulations and combined with continuum theory to pre- 
dict the macroscopic fracture energy 0. Here we focus 
on the microscopic stress distribution and its relation to 

Smax ■ 



A. Disentanglement versus chain scission 

The craze can fail by two different mechanisms that 
depend on the chain length N: short chains can disen- 
tangle, while very long chains fail through chain scission 
Both limiting behaviors and the crossover between 
them can be addressed through our simulations. As can 
be seen in Fig. short chains of length N = 128 form 
crazes that grow at the constant plateau stress S, but the 
stress drops monotonically to zero upon straining past A. 
For longer chains, the stress cr zz rises to a maximum value 
S'max that exceeds S. 

Values for S ma x were systematically obtained as a func- 
tion of normalized chain length N/N e from curves such 
as those shown in Fig. El Fig. 1221 summarizes the break- 
ing stresses for the craze fibrils normalized by the break- 



ing stress in the limit of very long chains Sqo. S ma x is 
zero for N < 2N e , since stable crazes do not form for 
such short chains. S m ax/Soo first rises roughly linearly 
with N/N e , and then saturates at unity for chain lengths 
longer than about 10iV e . The saturation coincides with 
the observation of significant amounts of chain scission. 
Interestingly, the data seems to collapse onto a single 
curve (solid line). 

Note that the maxima of the stress-strain curves in 
Fig. [7] are reached at strains of ~ 6 and ~ 10 for flexible 
and semifiexible chains, respectively. These values are 
close to v3A max , which implies that at the breaking point 
the chains are pulled fully taut between entanglement 
points. This was confirmed by direct analysis of the craze 
microstructure. 



1. Chain end relaxation 

In order to understand the crossover regime and the 
competition between the two failure mechanisms, it is 
useful to study the distribution of tension along a given 
chain. Fig.l23lshows the tension as a function of distance 
from the chain end at several stages of craze breaking. 
Since the chain ends are identical, symmetry was used 
to improve statistics. In the unstrained craze (lowest 
curves), both flexible and semifiexible systems exhibit a 
constant stress in the center of the chain, but a relax- 
ation toward the free ends. The characteristic length 
scales for this relaxation were extracted by fitting an ex- 
ponential decay to the transition region. The values of 
the decay lengths N^ nd — 21 and N^ d — 13 are compa- 
rable to the characteristic length of N e /3 for stretched 
segments, but are not universal. Stronger adhesive inter- 
actions were found to increase . Upon straining the 
craze, the tension in the center of the chains and the val- 
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FIG. 23: Distribution of tension along the (a) flexible and (b) 
semiflexible chains (N=512). The lowest curves corresponds 
to the unstressed craze and the highest curves show the ten- 
sion at the breaking point. Two intermediate stages are also 
shown. The characteristic length scales iV cn d describe the end 
relaxation and were obtained by fitting the indicated part of 
the curves to an exponential relaxation. 



ues of rise. At the breaking point (last curves), the 

end relaxation extends over a length scale comparable to 
the entanglement length. 

These results help to formulate a simple argument for 
the universal curve plotted in Fig. 1221 The average dis- 
tance of an entanglement point from the chain end is 
N/4. We assume that the probability of disentanglement 
decreases exponentially with distance from the chain end, 
as suggested by the tension relaxation curves. The char- 
acteristic length scale iVend a * the breaking point in these 
curves was on the order of N e . Because of the above, 
we expect this length scale to be the characteristic de- 
cay length for the probability of disentanglement, and 
thus postulate a disentanglement probability of the form 
exp[-(N - 2Ne)/4N e ]. Here, N was reduced by 2N e , 
since for this chain length the disentanglement probabil- 
ity is one and the chain is free on either side. The max- 
imum stress can now be written as the limiting value of 
Soo times the probability for non-disentanglement, which 
gives 



Fig. [53 shows that this curve agrees well with the data. 



B. Global tension distribution 

The parameter governing chain scission and thus the 
value of iSoo is the distribution of tension in the polymer 
craze. In a previous paper [lOj , we reported that this dis- 
tribution is characterized by an exponential tail at large 
tensile forces, in analogy to jammed systems such as gran- 
ular materials |39| . This distribution is shown in Fig. 1241 
for flexible and semiflexible chains of length N = 512 
and several strain states. The tensile (positive) part of 
the distribution is well fitted by exp[— // (/)], where 
only the positive tensions are included in the average (/) . 
Note that (/) is the same for flexible and semiflexible 
chains at the plateau. The distribution with the steepest 
slope (smallest (/)) corresponds to the fully developed 
craze. Additional curves with higher (/) correspond to 
stressed crazes at the same strain with respect to the un- 
strained craze. Note that the semiflexible and flexible 
crazes have the same values of (/) at each strain. This is 
related to the fact that the stress So in the fibrils is in- 
dependent of N e (see Section lyTTl . The last curve shows 
the tension distribution at the breaking point where <r zz 
is largest (see also Fig. [7J . The effect of straining the 
craze only changes (/) , and all curves could be collapsed 
after normalizing by (/). 

We note that in our simulations, no scission or dis- 
entanglement occurs during craze growth proper. The 
fraction of bonds that break at a given average tension 
is Ar«/» = f™exp[-f/{f)]df = exp[-/ c /(/)]. In our 
simulations f c — 100/lj as described in Section ITTT1 The 
onset of scission can be estimated using a simple scaling 
argument: Nb on d s Pbr((f)) ~ 1- From this we estimate 
an average value of the tension at breaking 



(/) = fc/\nN bonds . 



(20) 



J Soo = 1 - exp[-(A^ - 2Ne)/4N e 



(19) 



For Nb on ds — 262144, this implies that scission will occur 
when (/) ~ 8.0/lj, which was confirmed by direct in- 
spection of the chains at the corresponding strains. Such 
high tensions only occur when the craze is strained past 
the extension ratio. The largest value of (/) observed 
with the present model during craze growth was 4.5/lj 
and occurred at very low temperature T = O.Oluo/fcs 
and r c = 2.2a. 

The degree of chain scission in experimental crazes is 
still a matter of debate, but it appears likely that at least 
some chains do experience scission. The absence of scis- 
sion in the present study is most likely due to the low 
monomeric friction coefficient of the bead-spring model. 
As the above argument showed, a relatively modest in- 
crease in the average tension will quickly lead to appre- 
ciable scission. More realistic polymer models should be 
able to capture this effect. An increase in system size 
will likewise raise the number of broken bonds. For a 
typical value of (/) = 3/lj 7 one bond would break for 
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FIG. 24: Distribution of tension in crazes of size 262144 beads 
with N=512 for (a) flexible and (b) semiflexible chains. Strain 
states correspond to the ones shown in Fig. 1231 The values for 
(/> are 2.9,4.9,7.6, and 11.8/lj, respectively. / LJ = 2Au /a 3 
is the breaking force of the LJ interaction. The straight lines 
correspond to exp[-//(/>]/(/). 



every 10 15 bonds. Note that the exponential tension dis- 
tribution leads to a logarithmic size dependence Eq. I|20() 
and allows for sequential bond breaking. The fibrils are 
thus much weaker than implied by the common simple 
assumption that all bonds carry the same tension and 
break when (/) = f c . 



VIII. SUMMARY AND CONCLUSIONS 

This paper presented molecular dynamics simulations 
of craze nucleation, widening and breakdown. Initial fail- 
ure of the LJ polymer glass occurred through shear in bi- 
axial loading. Only when all three principal stresses were 
tensile did cavitation and craze formation occur. How- 
ever, once past the nucleation phase, plane stress condi- 
tions are sufficient for continuing craze growth. Cavita- 
tional failure could be fitted to a cavitation criterion of 
the form T% ev = Tg + a c p. 

Craze widening proceeds in the simulations by a clearly 
identifiable fibril drawing process. This interpretation is 
also well-supported by experiments. The resulting craze 
microstructure is compellingly similar to TEM images of 
experimental crazes, and the length scales quantified by 



(D) and (Do) are within experimental limits. The simu- 
lations clearly establish the connection between extension 
ratio and entanglement length. In the glass, disentangle- 
ment is prevented and the entanglements act like chem- 
ical crosslinks. A microscopic analysis of the length of 
stretched chain segments has shown that unlike the case 
of a regular mesh, only a few segments are fully expanded 
to the entanglement length, and the average extension is 
only N e /3. The factor 1/3 arises from averaging over all 
angles that a given segment can form with the stretching 
direction. 

Another salient finding of this study is the exponential 
distribution of tension in the craze. The presence of large 
stress fluctuations makes chain scission much more likely 
than e. g. a Gaussian distribution or uniform loading. 
Since force distributions of this kind arc also often seen 
in conventional "jammed" systems such as foams, colloids 
and granular media, we have suggested |lfj| that a craze 
can be viewed as a system that jams under tension. 

The highly nonequilibrium nature of the force distri- 
bution, and the strong concentration of stress in the co- 
valent backbone bonds, formed the basis for our critique 
of the conventional capillary model of craze widening. 
The polymer glass is not a viscous fluid in the active 
zone, and the hydrodynamic description does not apply. 
The picture suggested by our simulations is that craz- 
ing is a form of localized shear deformation, but with a 
much greater mobilization of material than in standard 
shear yielding. The very similar rate and temperature 
dependence of S is another indication of the close rela- 
tion between the processes. Based on trends observed in 
the simulations, we have suggested that the local stress 
in the fibrils So = SX is independent of the entanglement 
length. Sq varies with temperature and strength of ad- 
hesive interaction in a manner very similar to the yield 
stresses for shear and cavitation. Establishing a precise 
connection between these characteristic stresses should 
be a most interesting direction for future work. 

A detailed analysis of the microstructure of crazes was 
also presented. The calculated structure factor is similar 
to measured scattering intensities. As in these experi- 
ments, a Porod analysis was used to extract a measure 
of the mean fibril diameter (D) from the structure fac- 
tor. While the extension ratio depends only on N e , the 
mean fibril diameter depends on many factors. The value 
of (D) increases with increasing T and with increasing 
strength of the van der Waals interactions. Chain stiff- 
ness has less effect, although (D) is larger for flexible 
chains than semiflexible chains at high temperatures. 

The distribution of fibril diameters was determined 
from the real space structure of the crazes. The average 
fibril diameter from this method was always smaller than 
that determined from the structure factor. The distribu- 
tion was also very wide with a variance that exceeded the 
mean and a tail extending to many times the mean. The 
radial distribution function for the fibrils shows almost 
no correlation, merely an exclusion minimum near the 
origin. Fibrils merge and split with each other directly, 
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rather than being joined by smaller cross-tie fibrils. 

The simulations described here capture the generic fea- 
tures of experiments on many different polymers and pro- 
vide previously inaccessible information about the dy- 
namics and microstructure. However, they are unable to 
address quantitative behavior of specific polymers. Fu- 
ture studies with chemically realistic potentials will be of 
great value, but require orders of magnitude more com- 
putational effort. 
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